#### In Text: Age Analysis ####

# Libraries
# library(tidyverse)
# library(rio)
# library(here)
# library(srvyr)
# library(survey)

# # data_pnas = import(here("Data","data_pnas.rds"))

# This computes expected support for all our variables, then subsets
# to just ignoring court decisions for 25/65 year olds

age_out = data_pnas |>
  mutate(age = year(date) - birthyr) |> 
  select(age, uid, weight, violence1re:violence6re, norm_judgesre:norm_loyaltyre) |> 
  pivot_longer(violence1re:norm_loyaltyre, names_to = "var", values_to = "val",
               values_transform = list(val = ~as.numeric(.x))) |>
  group_by(var) |> 
  nest() |> 
  mutate(preds = map(data, \(x){
    mod = svyglm(val ~ age, design = as_survey_design(x, ids = uid, weights = weight), family = "binomial")
    p = predict(mod,
                newdata = data.frame(age = min(x$age, na.rm = T):max(x$age, na.rm = T)),
                se.fit = T,
                type = "response") |> as.data.frame()
    data.frame(age = min(x$age, na.rm = T):max(x$age, na.rm = T),
               fit = p$response,
               se = p$SE)
  })) |> 
  select(var, preds) |> 
  unnest(preds) |> 
  filter(var == "norm_judgesre", age %in% c(25,65))

age_25 = filter(age_out, age == 25) |> pull(fit)
age_65 = filter(age_out, age == 65) |> pull(fit)

cat(round(age_25*100,1), "% of 25-year olds and ",
    round(age_65*100,1), "% of 65-year olds support ignoring out-party court decisions",
    sep = "")